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The  Lagrange- Galerkin  method  for  the  two-dimensional 
shallow  water  equations  on  adaptive  grids 


Francis  X.  Giraldo*’1 


Naval  Research  Laboratory ,  Monterey,  CA ,  U.S.A. 


SUMMARY 

The  weak  Lagrange -Galerkin  finite  element  method  for  the  two-dimensional  shallow  water  equations  on 
adaptive  unstructured  grids  is  presented.  The  equations  are  written  in  conservation  form  and  the  domains 
are  discretized  using  triangular  elements.  Lagrangian  methods  integrate  the  governing  equations  along 
the  characteristic  curves,  thus  being  well  suited  for  resolving  the  non-linearities  introduced  by  the 
advection  operator  of  the  fluid  dynamics  equations.  An  additional  fortuitous  consequence  of  using 
Lagrangian  methods  is  that  the  resulting  spatial  operator  is  self-adjoint,  thereby  justifying  the  use  of  a 
Galerkin  formulation;  this  formulation  has  been  proven  to  be  optimal  for  such  differential  operators.  The 
weak  Lagrange -Galerkin  method  automatically  takes  into  account  the  dilation  of  the  control  volume, 
thereby  resulting  in  a  conservative  scheme.  The  use  of  linear  triangular  elements  permits  the  construction 
of  accurate  (by  virtue  of  the  second-order  spatial  and  temporal  accuracies  of  the  scheme)  and  efficient  (by 
virtue  of  the  less  stringent  Courant-Friedrich-Lewy  (CFL)  condition  of  Lagrangian  methods)  schemes 
on  adaptive  unstructured  triangular  grids.  Lagrangian  methods  are  natural  candidates  for  use  with 
adaptive  unstructured  grids  because  the  resolution  of  the  grid  can  be  increased  without  having  to  decrease 
the  time  step  in  order  to  satisfy  stability.  An  advancing  front  adaptive  unstructured  triangular  mesh 
generator  is  presented.  The  highlight  of  this  algorithm  is  that  the  weak  Lagrange- Galerkin  method  is 
used  to  project  the  conservation  variables  from  the  old  mesh  onto  the  newly  adapted  mesh.  In  addition, 
two  new  schemes  for  computing  the  characteristic  curves  are  presented:  a  composite  mid-point  rule  and 
a  general  family  of  Runge-Kutta  schemes.  Results  for  the  two-dimensional  advection  equation  with  and 
without  time-dependent  velocity  fields  are  illustrated  to  confirm  the  accuracy  of  the  particle  trajectories. 
Results  for  the  two-dimensional  shallow  water  equations  on  a  non-linear  soliton  wave  are  presented  to 
illustrate  the  power  and  flexibility  of  this  strategy.  Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 
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1.  INTRODUCTION 

The  advection  terms  in  the  governing  equations  of  fluid  motion  present  formidable  challenges 
to  many  spatial  discretization  methods,  including  Galerkin  methods.  These  terms  prevent  the 
operator  from  being  self-adjoint  and  as  a  result,  Galerkin  methods  are  no  longer  optimal  for 
the  spatial  discretization.  Researchers  have  tried  circumventing  this  problem  by  using  high- 
order  Eulerian  methods  and  characteristic -based  Lagrangian  methods.  In  this  section  a  brief 
overview  of  some  of  the  more  interesting  methods  is  given,  identifying  their  weaknesses  and 
assessing  their  strengths. 

A  limited  number  of  new  Eulerian  methods  have  been  introduced  in  recent  years  for  solving 
hyperbolic  partial  differential  equations  (PDEs)  (such  as  the  shallow  water  equations),  perhaps 
the  best  being  the  spectral  element  method  [1-3].  This  method  offers  high  accuracy  solutions 
and  spectral  convergence  (provided  the  solution  is  smooth)  but  at  the  price  of  having  to  use 
small  time  steps  (due  to  the  explicit  solvers  typically  used)  and  structured  grids  (due  to  the 
restriction  that  the  elements  be  quadrilaterals).  Recently,  Giraldo  [4]  showed  how  to  increase 
the  time  step  by  combining  the  Lagrange-Galerkin  method  with  the  spectral  element  method, 
while  Taylor  and  Wingate  [5]  have  been  working  on  developing  a  triangular  spectral  element 
method. 

In  contrast,  there  exists  a  multitude  of  characteristic-based  methods.  Eulerian -Lagrangian 
methods  combine  the  standard  finite  element  procedure  with  the  method  of  characteristics.  By 
integrating  the  equation  in  time  along  the  characteristics  the  resulting  spatial  operator  becomes 
self-adjoint,  which  then  justifies  the  use  of  Galerkin  spatial  discretizations.  These  methods  can 
be  classified  into  the  following  three  groups:  the  direct  (or  strong)  Lagrange-Galerkin  method, 
the  weak  Lagrange-Galerkin  method,  and  the  Eulerian-Lagrangian  localized  adjoint  method. 
In  the  strong  Lagrange-Galerkin  method,  the  time  derivative  and  the  advection  terms  are 
combined  into  the  Lagrangian  derivative  and  the  resulting  operator  is  then  discretized  using 
the  standard  finite  element  method.  This  is  the  approach  used  by,  for  example,  Bercovier  and 
Pironneau  [6],  Bermejo  [7],  Douglas  and  Russell  [8],  and  Priestley  [9],  In  this  method  the  basis 
functions  are  the  typical  Lagrange  polynomials  used  in  the  standard  finite  element  method, 
which  are  only  dependent  on  the  spatial  co-ordinates.  The  success  of  this  method  hinges  on 
determining  the  values  at  the  feet  of  the  characteristics.  Related  to  this  method  are  the 
semi-Lagrangian  [10]  and  characteristic-Galerkin  [11]  methods  (see  also  the  Taylor- Galerkin 
method  [12]).  The  semi-Lagrangian  method  is  essentially  the  strong  Lagrange-Galerkin 
method  with  the  exception  that  the  spatial  discretization  is  achieved  through  finite  differencing; 
this  method  is  quite  ubiquitous  in  the  meteorology  community.  In  this  method,  the  values  of 
the  variables  at  the  feet  of  the  characteristics  are  obtained  through  interpolation.  The 
characteristic-Galerkin  method  avoids  interpolations  by  using  approximations  obtained  from 
Taylor  series  expansions.  Therefore,  the  computations  all  take  place  in  terms  of  the  Eulerian 
grid,  which  eliminates  the  difficulties  of  the  Lagrangian  grid  but  at  the  expense  of  having  to 
use  smaller  time  steps  due  to  the  stricter  stability  restrictions  that  govern  all  Eulerian  methods. 

In  contrast,  the  weak  Lagrange-Galerkin  method  uses  basis  functions  that  are  dependent  on 
both  space  and  time  by  forcing  the  basis  functions  to  disappear  along  the  adjoint  of  the 
transport  operator  (the  transport  operator  referring  to  a  homogeneous  conservation  law);  this 
results  in  the  basis  functions  vanishing  along  the  characteristics.  Using  the  conservation  form 
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of  the  governing  equations  simplifies  the  resulting  discretization  by  introducing  the  Reynolds’ 
transport  theorem  [13],  which  then  eliminates  the  boundary  terms  arising  from  the  integration 
by  parts  used  to  obtain  the  adjoint  of  the  transport  operator.  This  is  the  method  introduced 
by  Benque  et  al.  [14,15]  and  the  approach  used  in  this  paper. 

Finally,  the  Eulerian-Lagrangian  localized  adjoint  method  forces  the  basis  functions  to 
disappear  along  the  adjoint  of  the  entire  differential  operator  (not  just  the  transport  operator). 
Since  the  Reynolds’  transport  theorem  is  not  used  in  this  formulation,  all  of  the  boundary 
terms  remain,  which  makes  it  rather  easy  to  enforce  any  type  of  boundary  conditions.  This  is 
the  method  introduced  by  Celia  et  al.  [16]  and  used  in  References  [17-20].  Thus,  for 
applications  where  the  boundaries  play  an  important  role  it  would  seem  prudent  to  use  this 
method.  However,  since  the  final  objective  of  our  work  is  to  develop  a  shallow  water  model 
on  the  sphere,  where  there  are  no  boundary  conditions,  the  weak  Lagrange -Galerkin  method 
appears  to  be  the  perfect  candidate  for  such  a  model. 

Lagrange -Galerkin  methods  have  increased  in  popularity  in  the  last  10  years  because  they 
offer  increased  accuracy  and  efficiency  by  virtue  of  their  independence  on  the  Courant- 
Friedrich-Lewy  (CFL)  condition;  however,  almost  all  research  has  involved  the  strong  method 
as  opposed  to  the  weak  method  presented  here.  Currently  the  only  applications  (of  which 
this  author  is  aware)  using  the  weak  method  or  variations  thereof  are  the  following:  two- 
dimensional  advection  [14,21],  two-dimensional  shallow  water  on  the  plane  [22],  two- 
dimensional  Navier-Stokes  [15],  and  two-dimensional  shallow  water  on  the  sphere  [23], 

There  have  been  few  attempts  at  using  the  weak  Lagrange -Galerkin  method  for  the  shallow 
water  equations.  In  Reference  [22],  the  weak  Lagrange-Galerkin  method  of  Benque  [15]  is 
presented  for  estuary  flows  on  structured  quadrilateral  grids.  In  their  paper,  the  interpolations 
for  the  departure  point  values  at  the  vertices  of  the  elements  are  bilinear  and  the  variables 
inside  their  respective  Lagrangian  element  are  then  assumed  continuous  and  linearly  varying 
throughout  the  element.  Assuming  a  linear  variation  within  the  Lagrangian  element  results  in 
a  highly  stable  scheme  but  a  strongly  diffusive  one  as  well.  In  our  paper,  triangular  grids  are 
used  and  the  integration  can  either  be  done  piecewise  exactly  or  numerically.  Neither  approach 
assumes  a  linear  variation  within  the  element  but  rather  takes  into  consideration  the  gradient 
across  the  different  Eulerian  element  that  the  Lagrangian  element  spans. 

In  Reference  [23]  an  exactly  integrated  weak  Lagrange-Galerkin  method  on  the  spherical 
shallow  water  equations  is  presented.  The  basis  functions  used  are  not  the  natural  co-ordinates 
derived  for  linear  triangles.  In  our  paper,  the  natural  co-ordinates  are  used.  This  means  that 
we  can  not  only  derive  explicit  basis  function  formulas  that  can  be  used  for  computing 
derivatives  or  searching  for  inclusion  of  departure  points,  but  also  they  can  be  used  to 
construct  exact  integration  formulas  for  all  of  the  finite  element  integrals  including  those 
occurring  at  the  Lagrangian  elements.  In  addition,  the  use  of  the  explicit  natural  co-ordinates 
has  ramifications  for  higher-dimensional  models.  All  of  the  algorithms  presented  in  this  paper 
are  directly  generalizable  to  a  two-dimensional  spherical  or  three-dimensional  Cartesian 
shallow  water  model.  This  has  been  the  focus  of  our  previous  [21]  and  future  work.  Therefore, 
the  idea  of  using  the  natural  co-ordinates  on  the  triangle  in  two-dimensional  space  is  directly 
applicable  to  a  spherical  model.  For  a  full  three-dimensional  model,  rather  than  using  the 
linear  triangular  natural  co-ordinates,  the  linear  tetrahedral  natural  co-ordinates  are  used  with 
little  change  to  the  algorithms  presented  in  this  paper. 


Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Fluids  2000;  33:  789-832 


792 


F.  X.  GIRALDO 


The  method  presented  in  Reference  [22]  uses  quadrilaterals  instead  of  triangles,  which  are 
used  in  this  paper.  There  are  many  advantages  in  using  triangles  over  quadrilaterals.  Since  the 
triangle  is  the  two-simplex,  any  two-dimensional  domain  can  be  discretized  with  triangles  no 
matter  how  multiply-connected  (i.e.,  non-convex)  it  is.  In  addition,  it  is  irrelevant  to  the 
scheme  whether  the  flow  has  strong  shears  or  is  highly  rotational.  If  quadrilaterals  are  used, 
the  Lagrangian  elements  may  become  highly  skewed  or  completely  distorted,  whereas  this  is  of 
no  concern  with  triangles.  Using  quadrilateral  elements  also  requires  a  non-linear  iterative 
solver  in  order  to  compute  the  local  element  co-ordinates  from  the  global  co-ordinates 
obtained  from  the  trajectory  equation.  This  is  required  because  the  interpolations  within  an 
element  must  be  carried  out  using  the  local  co-ordinates  because  the  local  basis  functions  are 
defined  only  in  terms  of  these  co-ordinates  and  not  in  terms  of  the  global  co-ordinates.  On 
triangles,  because  the  basis  functions  are  defined  in  terms  of  the  global  co-ordinate,  this 
procedure  is  not  required.  Furthermore,  in  this  paper  we  show  how  to  apply  the  weak  method 
with  adaptive  unstructured  grids,  which  can  only  be  accomplished  with  triangles.  The  idea  of 
using  adaptive  unstructured  grids  with  Lagrange -Galerkin  methods  is  attractive  because  the 
resolution  of  the  grid  can  be  increased  without  the  need  to  decrease  the  time  step.  In  Eulerian 
methods,  the  time  step  must  be  decreased  in  order  to  satisfy  the  CFL  condition;  violating  the 
CFL  condition  results  in  instabilities  and  inaccuracies.  In  Lagrangian  methods,  on  the  other 
hand,  the  CFL  condition  is  less  stringent  and  typically  the  time  step  need  not  be  decreased. 
Numerical  tests  are  presented  that  show  the  effects  of  increasing  the  time  step  on  the  accuracy 
of  the  scheme. 


2.  SHALLOW  WATER  EQUATIONS 
The  two-dimensional  planar  shallow  water  equations  in  conservation  form  are 
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but  note  that  if  we  move  the  pressure  terms  to  the  right-hand  side,  we  get 
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This  system  can  now  be  written  more  compactly  as 
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where 
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,  and  S(f)  = 
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The  equations  are  solved  for  the  three  conservation  variables  tp,  q>u,  and  cpv.  Clearly,  we  could 
also  include  other  forcing  functions  within  S(f).  Before  proceeding  to  the  discretization  of  the 
weak  Lagrange -Galerkin  method,  let  us  first  look  at  the  discretization  obtained  by  an 
Eulerian  method.  This  will  serve  both  for  contrast  and  also  because  we  require  one  Eulerian 
time  step  prior  to  using  the  weak  Lagrange -Galerkin  method  because  it  requires  two  known 
time  steps. 


3.  EULER-GALERKIN  METHOD 


Beginning  with  Equation  (2)  we  can  define  an  Eulerian  finite  element  method  by  taking  the 
weak  form 


df 

at 


+  V  •  (fu)  —  S(f) 


dQ  =  0 


and  integrating  by  parts  (Green’s  theorem)  such  that 
V  •  (^fu)  =  i//\ •  (fu)  +  (fu)  •  \i// 
to  arrive  at 


P  -f 

^^dQ+l  n-u^fdT—  I  V^'(fu)dQ  — 
Jr  Jn 


S(f)  dQ  =  0 


Jn  St 

In  the  case  of  no-flux  boundary  conditions,  the  second  integral  vanishes.  In  other  words 
n-ui^f  dF  =  0 
and  we  are  left  with 


i//C—dQ.= 
a  St 


\\jj  •  (fu)  dQ  + 


S(f)  dQ 


The  resulting  system  of  ordinary  differential  equations  (ODEs)  can  now  be  written 
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dU 
d  t 


H(  U) 


which,  after  integrating  in  time  by  a  general  family  of  Range  Kutta  schemes,  yields 
U*  + 1  =  U”  +  AtpH(\J)k~ 1 
where 


P  =  -r, — TTT  ’  and  H(U)°  =  H(Ur 

M  —  k  +  1 

This  scheme  is  required  for  the  first  time  step  because,  as  you  will  see,  the  weak  Lagrange - 
Galerkin  method  although  a  two-time  level  scheme,  requires  two  known  times  in  order  to  be 
able  to  compute  the  particle  trajectories  accurately. 


4.  WEAK  LAGRANGE  GALERKIN  METHOD 


4.1.  Spatial  discretization 

The  weak  form  of  Equation  (2)  is 


* 


<3f 

—  +  V  •  (fu)  —  S(f) 
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dQ  =  0 


and  using  the  calculus  identities 


Jt^t>  =  ^Jt  +  f^di  and  V‘(^fu)  =  ^V‘(fu)  +  (fu)‘V^ 


and  integrating  by  parts  results  in 
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The  first  bracketed  term  of  Equation  (4),  using  Reynolds’  transport  theorem  [13],  can  be 
written  as 


d^ 
d  t 


(<H)  dQ 


8t 


m  +  \-Wfa) 


dQ 


and  the  second  bracketed  term  is  actually  the  characteristic  equation 


(5) 
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di// 
d  t 


dt 


+  u-V^  =  0 


(6) 


where 


(7) 


is  used  to  predict  the  particle  trajectories  along  which  the  basis  functions  vanish.  The 
characteristic  equation  is  equal  to  zero  because  we  are  constraining  the  basis  functions  to  be 
constant  along  the  particle  trajectories.  Therefore,  in  essence,  the  basis  functions  ^  are 
functions  of  both  space  and  time.  This  arises  from  only  considering  changes  of  the  conserva¬ 
tion  variables  along  the  trajectory  lines  (see  References  [15,22]).  In  addition,  in  Appendix  A, 
we  have  included  an  example  to  show  that  the  basis  functions  satisfy  the  characteristic 
equation. 

Substituting  Equations  (5)  and  (6)  into  Equation  (4)  yields  the  weak  Lagrange-  Galcrkin 
system 

^  £  m  dQ  =  £  ^S(f)  dQ  (8) 

Note  that  the  advection  operator  has  disappeared  from  the  equations;  however,  the  correct 
particle  trajectories  are  accounted  for  by  the  trajectory  equation  (7).  In  addition,  consider  that 
the  divergence  of  the  velocity  has  also  disappeared  or  rather  has  been  absorbed  by  Reynolds’ 
transport  theorem. 


Remark  1 

The  implicit  absence  of  the  divergence  term  is  deceiving  because  this  term,  which  controls  the 
dilation  (contraction  or  expansion)  of  a  fluid  volume,  is  accounted  for  by  virtue  of  the 
trajectory  equation.  Therefore,  if  the  fluid  volume  dilates  then  it  is  exactly  satisfied  by  the 
dilation  of  the  Lagrangian  element  through  the  projection  of  the  element  grid  point  trajecto¬ 
ries.  This  key  point  is  the  difference  between  the  weak  and  strong  Lagrange -Galerkin  methods 
(see  Reference  [21]  for  details)  and  is  the  reason  for  which  the  weak  method  conserves  better. 

Figure  1  shows  a  schematic  of  the  dilation  of  a  fluid  volume.  The  grid  points  £j,  E2,  and  77, 
represent  the  vertices  of  the  Eulerian  element  at  the  arrival  time  tn  +  l.  The  trajectory  equation 
is  then  solved  in  order  to  determine  the  particle  trajectories  of  these  three  grid  points.  Let  Lu 
L2,  and  7.,  represent  the  departure  points  of  the  Eulerian  grid  points.  In  other  words,  the 
points  L  correspond  to  where  the  points  E  resided  at  time  7".  Note  that  the  areas  of  the 
triangles  E  and  L  are  not  necessarily  equal,  and  this  is  where  the  dilation  of  the  control  volume 
is  accounted  for.  In  other  words,  since  the  trajectory  equation  contributes  the  advection 
portion  of  the  equations  this  operator  does  not  need  to  be  included  in  the  equations.  In 
addition,  since  the  change  in  areas  of  the  control  volumes  obtained  by  the  trajectory  equation 
contributes  the  divergence  normally  given  by  the  velocity,  then  this  term  also  does  not  need  to 
appear  in  the  equations. 
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Figure  1.  The  dilation  of  the  Eulerian  element  ( E )  area  to  the  Lagrangian  element  (L)  area  traced  by  the 
particle  trajectories.  The  numbers  1-6  refer  to  the  six  Eulerian  elements  surrounding  the  Lagrangian 

element. 


4.2.  Time  discretization 

Integrating  Equation  (8)  in  time  by  the  0  algorithm  yields 


(^f)dQ'!  +  1  = 


n"+ 1 


(ij/f)  dQ'1  +  At 


i//S(f)  dQ"  + 1  +  (1  -  0)  i//S(f)  dQ" 


n«  + 1 


(9) 


which  represents  a  two-time  level  scheme  and  gives  the  explicit  rectangle  rule,  the  trapezoid 
rule,  and  the  implicit  rectangle  rule  for  0  =  0,  5,  and  1  respectively.  In  this  paper,  0  =  k  is  used 
throughout  because  it  yields  a  second-order  accurate  scheme  and  is  unconditionally  stable  with 
respect  to  advection.  (Lagrangian  methods  are  not  unconditionally  stable  for  advection  with  a 
forcing  function,  though;  in  this  case,  the  Courant  number  must  be  chosen  in  order  to  satisfy 
ODE  stability  conditions,  which  are  much  less  restrictive  than  their  PDE  counterparts.)  The 
other  two  schemes  are  both  first-order  and  the  0  =  0  is  not  very  stable  while  the  9  =  1  is  too 
diffusive  (see  Reference  [24]  for  details). 
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The  trajectory  equation  (7)  is  required  for  closure.  Clearly,  Equations  (9)  and  (7)  define  a 
two-time  level  scheme  requiring  the  variables  at  times  tn  and  t"  '.  However,  in  the  following 
section,  it  is  shown  that  in  order  to  solve  the  trajectory  equation  (7)  accurately,  the  velocity 
field  at  previous  time  levels  is  required.  The  accurate  solution  of  the  trajectory  equation  is 
perhaps  the  most  important  part  of  Lagrangian  methods  and  is  in  fact  often  overlooked.  In  the 
next  sections,  we  compare  two  new  methods  for  calculating  the  particle  trajectories. 

4.2.1.  Composite  mid-point  rule.  In  most  Lagrangian  methods,  the  particle  trajectories  are 
computed  by  the  mid-point  rule 

x»  +  i  _x«  =  A/u(x"  +  1/2) 

which  requires  the  following  extrapolation  of  the  velocity  field 

u»  +  l/2=  3U»_  1  u»-l 
2  2 


In  addition,  since  we  need  to  know  x'!+1/2,  we  have  to  assume  that  the  trajectory  from 
x"  1  — >  x”  is  linear;  in  other  words 


+  1/2  __ 


vn  +  1 


+  X" 


2 


(10) 


which  then  yields  the  iterative  equation 

x»+l/2  =  x»+l  _^u(x«+l/2) 

We  begin  with  the  arrival  point  x'!+1  as  the  first  approximation  x"  +  1/2.  Continuing  in  this 
fashion  results  in  convergence  in  less  than  five  iterations.  This  scheme  is  second-order  accurate 
in  space  and  time  and  is  used  quite  extensively  in  the  semi-Lagrangian  community.  (The 
semi-Lagrangian  method  when  used  in  conjunction  with  the  finite  element  is  also  known  as  the 
strong  Lagrange -Galerkin  method.) 

However,  the  weakness  of  this  scheme  is  that  upon  obtaining  the  correct  mid-point 
trajectories,  the  departure  points  are  then  computed  assuming  a  linear  variation  in  the 
trajectories  as  shown  in  Equation  (10).  One  simple  way  of  improving  this  scheme  is  to 
construct  a  composite  mid-point  rule.  The  idea  is  to  sub-divide  the  trajectory  into  multiple 
trajectories.  Assuming  we  have  M  sub-divisions  of  the  trajectory,  we  begin  from  time  tn+ 1  and 
construct  the  following  scheme: 

x«  +  1  —  p  x»  4-  1  —  oc  _ u/  x«  -f-  1  — 

2 


where 
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i- 1  „  '-(1/2)  i 

a  = - ,  /)  = - ,  y  =  — 

M  '  M  r  M 


and  the  mid-point  departure  points  are 


+  1  —  y  _  ' 


+  1  —  a  _|_  +  1  —  /? 


for  i  =  1, .  . . ,  M 
now  given  by 


This  approach  defines  M  series  of  second-order  accurate  steps  in  both  space  and  time.  Note 
that  the  velocity  field  must  be  extrapolated  in  order  to  get  its  value  at  time  tn  + 1  -  p.  This  is 
achieved  by  the  following  second-order  extrapolation: 

u(?»+ '  -  /?)  _  (2  _  fj)u(  t")  -  (1  —  /])u(tn~l)  (11) 

which  then  requires  an  interpolation  to  get  the  values  u(x"  1  /;). 

4.2.2.  Runge-Kutta  scheme.  A  higher-order  approximation  can  be  obtained  by  applying  the 
Runge-Kutta  scheme  to  the  trajectory  equation  (7),  yielding 

xfc+1  =  x«  +  1  _  At/iu(x)k~  1 

where 


P 


1 

M  —  k  +  1  ’ 


k  =  1, .  . . ,  M  and  u(x) 


u(x)'' +  1  for  k  =  1 

u(x)"+l  for  k>  1 


This  requires  knowing  u(x"  1  which  can  be  extrapolated  from  the  velocity  fields  at  the 
previous  time  steps.  We  could  use  the  extrapolation  given  in  Equation  (11)  but  this  would 
formally  give  only  second-order  approximations  for  the  velocity.  Therefore,  the  most  consis¬ 
tent  approach  is  to  use  an  extrapolation  of  equal  order  to  that  of  the  Runge-Kutta  scheme. 
At  the  first  Lagrange -Galerkin  step  (tn)  we  already  know  the  velocity  at  times  /"  1  and  tn~2 
(because  the  first  time  step  tn~  1  is  carried  out  with  the  Eulerian  method  and  at  t"  2  the  intial 
conditions  are  used).  Therefore,  at  this  time  step  we  can  only  use  a  maximum  of  a 
second-order  scheme;  for  this  reason  we  use  a  second-order  extrapolation  in  conjunction  with 
a  second-order  Runge-Kutta  scheme  (M  =  2).  However,  at  the  next  time  step  (t'l  +  r)  we  know 
the  velocity  at  tn,  tn  ',  and  /"  2.  We  can  do  this  one  more  time  until  we  know  the  velocity 
field  at  sufficient  time  levels  in  order  to  obtain  a  fourth-order  scheme.  The  extrapolations  of 
the  velocity  field  at  time  /"  1  1  then  can  be  obtained  from  a  Taylor  series  expansion  giving 

a{tn+ 1  ~p)  =  (1  +  a  +  b  +  c)u(fM)  +  au(tn~  *)  +  Zr u(  r "  2 )  +  eu(?"~3)  (12) 
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where 


a  = 

8 

O 

II 

c 

=  0 

for 

M=  2 

a  = 

28 +  82 

b=  -\8 

-id2 

c 

=  0 

for 

M  =  3 

a  = 

38  +  §82  +  ±83 

b=  — 1<5 

-2d2  -id3 

c 

=  \5  +  \8z  +  \fi3 

for 

M=  4 

and 

8  =  \-P 

which  means  that  after  four  time  steps,  the  scheme  will  be  running  at  fourth-order  accuracy 
throughout  the  remainder  of  the  time  integration.  Note  that  this  scheme  does  not  require  a 
first  guess  starting  point  as  in  the  previous  scheme  because  this  is  not  an  iterative  scheme  but 
rather  a  multi-step  scheme.  Therefore,  we  begin  from  the  arrival  point,  as  in  the  previous 
scheme,  and  step  through  towards  the  mid-point.  The  advantage  of  this  scheme  over  the 
composite  mid-point  rule  is  that  the  Runge-Kutta  scheme  never  assumes  a  linear  trajectory, 
and  is  in  fact  inherently  non-linear.  The  superiority  of  this  scheme  over  the  composite 
mid-point  rule  becomes  evident  when  very  large  time  steps  are  taken.  Because  of  its  non-linear 
structure,  and  its  high-order  approximation  the  Runge-Kutta  scheme  is  better  equipped  to 
compute  more  accurate  trajectories  than  the  mid-point  rule  especially  when  large  changes  of 
the  velocity  field  with  respect  to  time  occur. 


4.3.  Element  equations 

Returning  to  Equation  (9)  and  substituting  in  the  conservation  variables  from  Equation  (3) 
gives 


/* 

V 

n  +  1 

r 

(p 

n 

* 

(pu 

dQ"  +  1=  >)/ 

cpu 

dQ" 

Jn«  +  i 

(pv 

Jn« 

(pv 

r 

0 

+  A  tO 

—  (p(d(p/dx)  +  ftpv 

JQP  +  1 

-  (p(8(p/dy)-f<pu 

n  +  1 


dQ' 


\n  +  1 


0 

f  ^ 

—  (p(8cp/8x)  +  ftpv 

Jq« 

_  -  <p(.d<p/dy ) 

dQ" 


(13) 


which  is  the  equation  to  be  solved  within  each  element.  The  mass  and  momentum  equations 
can  be  decoupled  as  follows: 
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the  mass  equation  is  given  as 


n«+ 1 


\j/((p )  dQ"  + 1  =  lP(<p)  dQ" 


which  is  no  longer  a  function  of  velocity  except  through  the  trajectory  equation, 
the  momentum  equations  are 


i//(<pu)  dQ'I+  1  —  A td 
cin  +  l  Jci"  + 1 

r  r 

i/t(<pu )  dQ"  +  At  (l  -  6) 


n>‘ 


Cl" 


i//(fcpv)  dQ"  + 1 

\J/(  —  cp  — ^  +  fcpv  I  dQ"  +  A tO 

ox 


Cl"  +  1 


and 


i//(<pv)  dQ"  + 1  +  A  tO 


Cl"+  1 


i//(fcpu)  dQ 


n  +  1 


Q«+  1 


i//((pv)  dQ"  +  Ar(  1  —  9) 


i/'l  —  cp  —  fcpu\  dQ"  +  AtO 

oy  }  Jci"+ 1 


*1  +  1 


where  they  are  coupled  through  the  Coriolis  terms.  Since  the  mass  equation  can  be  solved  for 
first,  the  gravity  terms  in  the  momentum  equations  ((p(d<p/dx)  and  (p(d(p/dy))  at  time  tn+ 1  can 
be  moved  to  the  right-hand  side  as  known  quantities.  These  derivatives  can  then  be  obtained 
by  differentiating  the  basis  functions  as  is  done  in  the  typical  finite  element  method.  However, 
these  same  gravity  terms  at  time  tn  need  to  be  determined  in  some  other  way  because  we  need 
these  values  at  the  quadrature  points.  Section  4.5  describes  how  these  derivatives  are  obtained. 
Approximating  the  conservation  variables  within  each  element  Q  by  the  relation 

re)=  i  t jfj  (Hi 

j=  1 

gives  the  following  matrix  form  for  the  mass: 


[MA  dQ"+  [{<p) 


1/7+1  _ 


Q«+  1 


M<?)]”dQ" 


and  for  the  two  momentum  equations 
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» 

'I'tl'j 

-  Atd\/ijlfjliJk 

dQ"  + 1 

( <PU)j 

Jn«  +  i 

_  +  A  t6^^ij4Jk 

'I'i'l'j 

>l/i((pu)  +  A/(l  -  6)1/// 
i/Ziicpv)  +  At(l  -6)1//, -I 


d<P  ,  r 

~(p  —  +f(pv 

OX 

dtp 

-v-r 

dy 


n 


clQ" 


+ 


Ja«  + 1 


-  A tOi/z^jcpj  —  cpk 

-  A tOij/^jcpj  —  cpk 


n  +  1 


dQK  +  1 


(15) 


which  represents  the  element  system.  In  Section  4.5,  we  describe  how  the  integrals  at  time  tn 
are  obtained. 

4.4.  Global  equations 

Adding  all  of  the  contributions  from  the  elements  to  each  of  the  global  grid  points  results  in 
the  following  linear  system  of  global  equations  for  the  mass 

Aij(p]+l=bf  (16) 

and  the  momentum  equations 


~Aj 

-  Bu 

((pu)j 

n  +  1 

b / 

L  Bej 

Aj  J 

_K_ 

where 


A,  ,=  |  ^.dQ"  +  1 


Y  IT  J  • 

Jn«  + 1 


blj  =  Ate 


n«+i 


1 1'i'/'ji/'kfk  dfi' 


i«  +  1 


(17) 


and  bf,  h ",  and  b1/  are  the  right-hand  side  vectors  for  the  mass  and  u  -v  momentum 
respectively. 

The  matrices  A  and  B  are  both  symmetric  positive  definite  and  so  the  mass  matrix  is  also 
symmetric  positive  definite  but  the  momentum  matrix  is  not.  The  mass  matrix  poses  no 
difficulty  and  can  be  solved  quite  rapidly  using  conjugate  gradient  methods.  However,  the 
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momentum  matrix  is  skewed  symmetric  and  for  this  type  of  matrix  there  is  no  one  best  solver. 
The  momentum  matrix  is  sparse,  however,  and  iterative  methods  may  be  the  best  methods  of 
solution.  The  good  news,  however,  is  that  although  we  began  with  a  system  of  non-linear 
equations,  discretizing  the  equations  by  the  weak  Lagrange -Galerkin  method  results  in  a 
system  of  decoupled  linear  equations. 

Because  the  momentum  matrix  is  not  symmetric  positive  definite  we  cannot  safely  utilize 
many  of  the  best  known  solvers,  such  as  the  conjugate  gradient  method.  There  are,  however, 
variants  of  these  methods  that  can  be  used,  such  as  the  bi-conjugate  gradient  method.  Because 
the  system  is  quite  sparse,  direct  methods  are  inefficient.  For  the  moment,  we  are  employing 
a  direct  lower-upper  (LU)  decomposition  method,  which,  while  not  the  most  efficient,  offers 
a  simple  way  of  inverting  the  coefficient  matrices. 

Another  possibility  for  solving  the  momentum  equations  more  efficiently  is  to  simplify  the 
equations  in  order  to  make  them  symmetric  positive  definite.  This  approach  is  introduced  in 
Reference  [23]  and  the  strategy  is  to  extrapolate  the  velocity  field  (using  the  times  t"  and  t" ~ ') 
to  get  an  approximation  at  tn+l.  This  now  allows  moving  the  B  matrices  in  Equation  (17)  to 
the  right-hand  side.  Each  conservation  variable  can  now  be  solved  explicitly.  Flowever,  doing 
this  greatly  restricts  the  maximum  time  step  that  can  be  used.  In  our  paper,  we  invert  the 
skewed  symmetric  system  with  a  direct  solver  but  are  currently  exploring  more  efficient 
iterative  methods,  such  as  the  multi-grid  and  bi-conjugate  gradient  methods. 

4.5.  Basic  functions 

The  conservation  variables  belong  to  the  following  spaces: 

(peHl(Cl),  ( (pu ,  (pv)eHl(Cl ) 
and  their  test  functions,  likewise,  are  defined  as 
ffeH\ Q), 

in  other  words,  they  belong  to  the  set  of  square  integrable  functions  whose  first  derivatives  are 
also  square  integrable  and  are  the  linear  triangular  basis  functions.  The  Sobolev  space  H o(Q) 
is  defined  as 


//o(Q)  =  {^e77o(Q)|^(r)  =0} 


where  T  denotes  the  boundary  of  the  domain  Q.  These  functions  are  written  as  follows: 


atx  +  bj’  +  ct 
det  AiJJc 


and  have  the  exact  integration  rule 


» 

Jq 


fl'I'j'Pl  d-Q  =  det  AiJJc 


<x\p\y\ 

(a  +  f  +  y  +  2)1 


(18) 


(19) 
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where 

det  AiJJc  =  (xj  -  x,.)  x  (xk  -  x,  ) 


and 


at  =  yj  -  yk ,  bt  =  xk  -  Xj,  c,  =  Xjyk  -  xkyj 

The  indices  i,j,k,  =  1, .  . . ,  3  are  cyclic  in  the  following  order: 

k 


>/  \ 
i  ->  j 

The  C°  basis  functions  given  in  Equation  (18)  are  also  called  natural  or  area  co-ordinates 
because  they  represent  the  area  ratios  between  the  triangle  Aijk  and  the  sub-triangles  formed 
by  any  point  p  in  the  interior  of  the  element  and  the  three  vertices  i,j,  and  k.  For  this  reason, 
using  these  types  of  functions  extend  to  higher  dimensions  because  we  can  use  similar  triangles 
on  the  sphere  as  in  References  [21,25],  or  we  can  use  the  volume  co-ordinates  that  are  the 
linear  tetrahedral  functions  in  R3. 

4.5.1.  Derivatives.  The  approximation  of  any  function  /  within  an  element  is  given  by  Equation 
(14)  and  its  derivative  in  the  ^-direction  follows  from  Equation  (14)  and  can  be  written  as 


8s  jh  8s  Jj 


(20) 


where  s  can  be  either  the  x  or  y  co-ordinate.  Recall  that  in  Equation  (15)  we  need  to  know  the 
values  of  the  gravity  terms  <p(d<p/ds )  both  at  times  tn+ 1  and  tn.  Clearly,  at  time  tn+  1  we  can 
use  Equation  (20)  because  these  integrals  occur  in  the  Eulerian  elements  (see  Figure  1). 
However,  at  time  t'\  integrals  now  occur  in  the  Lagrangian  elements  and  if  these  elements  span 
more  than  one  Eulerian  element,  then  using  (20)  would  result  in  a  constant  derivative 
approximation  across  elements,  which  is  clearly  not  consistent  with  C°  basis  functions. 
Therefore,  the  best  way  to  obtain  these  values  is  the  same  way  in  which  any  other  value  at  the 
quadrature  points  is  obtained;  i.e.,  by  using  the  vertex  (or  node)  values  of  the  Eulerian 
elements  to  interpolate  (using  the  local  linear  basis  functions)  the  values  at  the  quadrature 
points.  For  the  conservation  variables  <p,  (pit,  cpv,  these  values  are  known  at  all  the  grid  points; 
however,  the  derivatives  are  not.  One  way  of  obtaining  these  grid  point  derivatives  is  to 
consider  the  following:  the  derivative  within  an  element  can  be  written  as  done  previously 
using  the  function  values  themselves  at  the  grid  points  as  given  in  Equation  (20).  However,  if 
the  derivatives  at  the  grid  points  themselves  were  known,  then  we  could  also  write  the 
derivatives  within  an  element  as 
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8s 


% 

8s 


and  equating  Equations  (20)  and  (21)  in  a  weak  sense,  yields 


Jn 


'I'i'l'i  dQ 


Vj 

8s 


(21) 


(22) 


where  the  element  contributions  are  summed  in  order  to  form  a  global  system,  which  when 
inverted  yields  the  derivatives  at  the  grid  points.  This  derivative  matrix  is  symmetric  positive 
definite  and  can  be  inverted  easily;  however,  it  can  also  be  diagonalized  for  the  sake  of 
efficiency.  This  is  the  approach  used  in  our  paper.  In  Reference  [25]  this  derivative  computa¬ 
tion  approach  was  shown  to  be  second-order  accurate  and  better  than  a  centered  finite 
differencing  approximation. 

4.5.2.  Integrals.  At  the  Eulerian  grid  points  (those  at  time  tn+ ')  the  integration  is  done  exactly 
using  Equation  (19)  but  it  still  remains  to  be  decided  how  to  handle  the  Lagrangian  elements 
(those  at  time  t").  One  way  is  to  write  the  approximation  within  the  Lagrangian  element  and 
equate  the  values  using  the  Eulerian  element,  which  it  intersects.  Let  /  be  the  region  LrsE. 
Thus,  we  can  write  any  grid  point  shared  by  both  elements  as 

x  =  i  <a/x/=  i  ry 

j= 1  7=1 

and  expanding  these  relations  results  in  the  following  matrix  system: 


x[  x{  x{ 

>r 

'xftf  +  xfyZ  +  xftf 

y\  y{  y\ 

= 

yfi^f  +  y^  +  yUf 

l  l  l 

rK 

1 

which  after  inverting  yields 

o#f  +c#f+4 

*■' - - 

where 

a^xfni  +  yfa  b^xfri+yK1, 
£!=xl-x',  n\=yIj—yIk 

and 

det  A O*  =  (x/  -  xf)  x  (x(  -  xj) 


(23) 


c\  =  xf >/f  +  yUl  4  =  x/y(  -  xfy / 
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In  Figure  1  we  see  that  the  Lagrangian  element  Ll  2, 3  spans  six  Eulerian  elements.  Note  that 
element  1  and  the  Lagrangian  element  overlap  in  only  one  triangular  intersection  region  (dark 
triangle),  whereas  element  2  overlaps  with  the  Lagrangian  element  in  two  triangular  intersec¬ 
tion  regions  (light  triangles).  The  basis  functions  of  these  intersection  regions  are  written  as 
functions  of  the  surrounding  Eulerian  basis  functions.  Since  we  have  exact  integration  rules  for 
any  integral  involving  the  Eulerian  basis  functions,  we  now  also  have  exact  integration  rules 
for  the  intersection  regions.  This  exact  integration  approach  was  introduced  in  Reference  [9] 
and  is  quite  an  elegant  method.  Unfortunately,  in  order  to  integrate  the  elements  exactly  we 
must  first  decompose  the  Lagrangian  element  such  that  the  pieces  lie  exactly  within  Eulerian 
elements.  This  requires  an  unstructured  mesh  generation  strategy,  which  essentially  becomes  a 
very  tedious  exercise  in  intersection  and  inclusion  tests.  Therefore,  in  the  absence  of  efficient 
intersection  and  inclusion  tests  it  is  best  to  use  some  other  simpler  means. 

Another  approach  for  obtaining  the  integrals  at  time  t"  is  simply  to  interpolate  the  vertex  of 
the  Lagrangian  element  (Lu  L2,  L3)  and  then  construct  linear  basis  functions  within  this 
element.  The  obvious  shortcoming  of  this  approach  is  that  we  now  lose  the  variation  within 
the  element  by  assuming  a  linear  approximation.  Using  this  simpler  approach  would  render  the 
method  quite  inutile.  Therefore,  the  final  possibility  is  to  take  sufficient  sampling  (or 
quadrature)  points  in  order  to  capture  the  variation  within  the  element.  In  other  words,  we 
seek  the  minimum  number  of  quadrature  points  that  can  extract  the  maximum  amount  of 
information  from  all  six  surrounding  elements. 

Many  quadrature  rules  have  been  tested  and  the  quintic  (seven-point)  Gauss  quadrature  rule 
appears  to  be  the  minimum  that  should  be  used  in  order  to  get  results  comparable  with  the 
exact  integration  approach.  In  the  numerical  integration  approach,  the  vertices  of  the  Eulerian 
triangular  element  (tn+l)  are  mapped  along  the  characteristics  to  where  they  resided  at  the 
previous  time  step  ( tn ).  The  feet  of  these  characteristics  now  forms  the  Lagrangian  element 
Lx  2  3*  Within  this  Lagrangian  element  the  variables  at  the  quadrature  points  are  interpolated 
linearly  using  the  basis  functions  (area  co-ordinates)  of  the  Eulerian  element,  which  claims  the 
quadrature  points.  In  other  words,  the  quadrature  value  arising  from  the  Eulerian  element  1 
(Figure  1)  is  obtained  by  linearly  interpolating  the  vertex  values  of  element  1  at  the  quadrature 
point.  Similarly,  the  quadrature  point  lying  in  element  2  gets  its  value  from  the  interpolation 
of  the  values  from  element  2,  and  so  on.  In  this  fashion,  the  integration  takes  into  account  the 
variation  (higher  than  linear)  of  the  function  values  inside  the  Lagrangian  element.  For  this 
reason,  using  a  low-order  quadrature  rule  will  result  in  large  inaccuracies  while  higher  orders 
will  always  result  in  increased  accuracy  (until  a  limit  is  reached).  For  the  problems  studied  in 
this  paper,  the  minimum  number  of  sampling  points  appeared  to  be  around  seven,  which 
results  in  our  choice  of  a  quintic  quadrature  rule. 

Using  numerical  integration  results  in  almost  identical  results  to  those  of  exact  integration. 
The  only  potential  pitfalls  are  in  conservation  and  stability.  The  exact  integration  approach  is 
exactly  conservative  provided  that  a  direct  solver  is  used  to  invert  the  matrix.  This  is  rarely 
feasible,  and  when  iterative  methods  are  used  exact  conservation  can  no  longer  be  guaranteed. 
In  practice,  however,  the  lack  of  exact  conservation  is  insignificant  and  has  not  been  shown  to 
adversely  affect  the  results.  For  the  problems  studied  thus  far,  conservation  has  been  preserved 
even  for  the  numerical  integration  approach  with  an  iterative  solver. 
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The  question  of  stability  is  a  bit  more  elusive  to  pin  down.  Priestley  [9]  showed  that  using 
quadrature  rules  for  the  numerical  integration  of  the  Lagrangian  elements  results  in  instabili¬ 
ties,  in  theory.  However,  in  practice  these  instabilities  have  yet  to  present  themselves  either  in 
Priestley’s  [9]  or  our  current  work.  The  exact  integration  approach  does  not  suffer  the 
likelihood  of  such  instabilities,  however.  A  comparison  of  the  exact  and  numerical  integration 
approaches  are  discussed  in  Section  6. 


4.6.  Searching  algorithms 

The  crux  of  any  Lagrangian  scheme  is  the  accurate  calculation  of  the  particle  trajectories,  in 
other  words  the  correct  solution  of  Equation  (7).  Once  the  correct  trajectories  are  computed  it 
then  remains  to  interpolate  either  the  quadrature  points  in  the  numerical  integration  approach 
or  the  points  comprising  the  intersection  regions  in  the  exact  integration  approach.  For 
structured  or  quasi-structured  grids  ad  hoc  searching  algorithms  can  be  constructed  but  for 
unstructured  or  general  grids  quadtree  algorithms  are  the  best  choice. 

The  idea  behind  the  quadtree  algorithm  is  to  find  the  nearest  neighbor  in  the  background 
(Eulerian)  grid  to  the  departure  point  being  searched.  Let  quad_tree[  1:  ntree,  1:  7]  be  an 
integer  array  that  stores  this  quadtree.  This  array  stores  the  following  information: 

•  quad_tree[i,  1-4]  stores  the  four  children  of  this  quad. 

•  quad_tree[i,  5]  stores  the  position  of  this  quad  with  respect  of  its  parent. 

•  quad_tree[i,  6]  stores  the  location  of  its  parent. 

•  quad_tree[i,  7]  stores  the  number  of  nodes  contained  within  this  quad. 

However,  another  integer  array  is  required,  which  stores  the  list  of  elements  claiming  each 
grid  point.  Upon  finding  this  nearest  neighbor  we  then  search  through  the  list  of  elements  that 
claim  this  grid  point  and  check  for  inclusion  using  the  local  area  co-ordinates  (basis  functions) 
of  the  element.  If  we  let  the  departure  point  be  denoted  by  p,  then  if  the  following  conditions 
are  met: 


det  Ap j' k  >  0,  detAii)  fc>0,  and  det  dsi  jp  >  0 

then  the  point  p  is  said  to  be  contained  within  the  element  with  vertices  i,  j ,  k. 

There  are  usually  no  more  than  six  elements  claiming  each  node  even  for  distorted 
unstructured  grids.  For  highly  distorted  grids,  however,  the  departure  point  may  not  necessar¬ 
ily  lie  within  one  of  the  elements  claiming  the  nearest  neighbor.  A  good  example  where  this  can 
occur  is  in  non-convex  domains.  In  this  case,  during  the  sweep  through  the  element  list 
claiming  the  nearest  neighbor,  the  minimum  distance  between  the  departure  point  and  the 
element  nodes  is  stored.  The  element  node  yielding  the  minimum  distance  is  considered  to  be 
the  new  nearest  neighbor.  If  no  inclusion  is  found,  then  the  new  nearest  neighbor  is  used  and 
the  process  is  repeated.  Therefore,  in  the  worst  case  scenario,  only  two  nearest  neighbor  loops 
are  required.  This  can  have  adverse  affects  on  the  efficiency  of  the  scheme  if  this  case  should 
arise  often;  fortunately,  this  situation  has  not  yet  presented  itself. 
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5.  ADAPTIVE  UNSTRUCTURED  MESH  GENERATOR 


A  practical  condition  for  creating  optimal  grids  based  on  numerical  solutions  is  to  require  that 
the  error  be  equally  distributed  across  all  elements  fi  comprising  the  domain  /).  Let 

N 

cpi(x,y)=  X  ^ m(x,y)(p{xm,ym ) 

m  —  1 


be  the  interpolation  of  <p(x,  y)  in  Q,  where  N  is  the  number  of  interpolation  points.  The  error 
in  using  the  interpolant  can  be  written  via  a  Taylor  series  expansion  about  the  point  (a,  b)  as 


cp (x,  y)  -  q>j(x,  y )  =  £  <p (iJ\x,  y) (A  U)  -  b) 

i+j=k+  i  ?!  j'- 


-  X  X  <PiiJ\x,y) 

m  —  1  i  +  j  =  k  +  1 


(x  —  ay  (y  —  b)J 

j'- 


yjx,y) 


where  cp(iJ>  denotes  the  derivatives  of  cp  with  respect  to  x  and  y.  Taking  the  oo  norm  and  using 
norm  properties  to  simplify  gives 


\<p(x,  y)  —  (Pi(x,  y)  j|  oo  <  C  max  hk  +  1  X  \(p(iJ\x,  y)\  (24) 

i+j=k+ 1 

In  this  relation  k  is  the  order  of  the  interpolation  functions,  h  is  the  local  grid  spacing,  and  C 
is  a  constant.  A  value  of  k  =  1  corresponds  to  linear  triangles.  Therefore,  for  an  optimal  grid 
using  linear  triangular  finite  elements  we  require  that 

Cmax/i2  X  |^(,’-/)(x,  y)\  =  costant,  VQ  eD 

i+j  —  2 

The  choice  of  error  indicator  and  indicator  variable  are  very  much  problem-dependent.  The 
error  indicator  used  in  this  paper  uses  the  derivative  of  the  mass  variable.  The  strategy  behind 
this  mesh  generator  can  be  characterized  as  a  global  reconstruction  (or  remeshing),  as  opposed 
to  a  refinement  method  that  only  changes  the  grid  locally.  It  is  an  advancing  front  method  and 
is  termed  a  global  reconstruction  because  whenever  the  grid  must  be  altered  it  is  constructed 
beginning  with  the  boundaries  and  circling  towards  the  center  of  the  domain  (see  Reference 
[26]  for  further  details  on  this  mesh  generator). 

The  accurate  interpolation  of  values  from  one  grid  to  a  newly  refined  grid  is  not  at  all  simple 
since  both  grids  (the  old  and  newly  adapted  meshes)  may  be  completely  unstructured  and  the 
domain  is  not  necessarily  convex.  Therefore,  local  methods  must  be  used  whereas  global 
constructs  like  splines  cannot. 

For  steady  state  problems,  the  grid  points  on  the  new  mesh  can  be  interpolated  using  the 
basis  functions  and  elements  from  the  old  mesh.  This  is  possible  because  we  only  require  an 
initial  guess  for  the  problem  and  while  the  initial  guess  can  affect  the  convergence  of  the 
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solution,  it  does  not  affect  the  accuracy.  In  transient  problems,  on  the  other  hand,  the  new  grid 
point  values  should  not  only  be  as  accurate  as  possible  but  the  variables  should  also  be 
conserved.  By  using  the  idea  of  the  weak  Lagrange -Galerkin  method,  we  can  obtain  these 
conserving  integrals,  which  yield  the  values  at  the  new  grid  points  [9],  In  other  words,  we  solve 
the  following  system: 


dQnew<p“ew 


^olddQ0 


^old 


where  the  right-hand  side  is  known.  This  approach  globally  conserves  the  variable  anywhere 
from  one  to  three  orders  of  magnitude  better  than  a  linear  interpolator  using  the  local  basis 
functions.  In  fact,  the  change  in  the  conservation  variables  from  an  old  mesh  projected  onto 
a  newly  adapted  mesh  is  around  the  order  of  1  x  10  " 8  percent,  which  is  obviously  conservative 
to  within  machine  precision. 


6.  NUMERICAL  EXPERIMENTS 


For  the  numerical  experiments,  the  following  terms  are  used  to  compare  the  performance  of 
the  schemes:  the  L2  error  norm 


(/exact  -ff  dU 
/exact  d “  “ 

Jn 


where  /  represents  any  of  the  conservation  variables  in  Equations  (3),  the  trajectory  norm  [25] 


(xD  -  x“act)2  dQ 
r=  jn - 

(xA  -  xgact)2  dQ 
Jn 

where  subscripts  A  and  D  denote  the  arrival  and  departure  points,  and  the  following  two 
additional  measures 


M  = 


q>  dQ 


exact  *JQ 


E  = 


cp(u2  +  v2)  +  cp 2  dQ 
Jo 

/» 

^exact(^  exact  "h  ^  exact)  T  exact 

Jq 


The  L2  error  norm  compares  the  root-mean-square  percent  error  of  the  numerical  and  exact 
solutions,  T  measures  the  accuracy  of  the  trajectories,  M  measures  the  conservation  property 
of  the  mass,  and  E  measures  the  conservation  of  the  total  available  energy.  The  ideal  scheme 
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should  yield  L2  error  and  trajectory  norms  of  zero,  and  mass  and  energy  measures  of  one.  In 
addition,  the  following  ratio  is  used  throughout  the  following  sections: 

At 


The  variable  represents  the  ratio  between  the  Lagrange -Galerkin  time  step  and  the  maximum 
allowable  time  step  for  the  Eulerian  fourth-order  Runge-Kutta  presented  in  Section  3.  Since 
the  fourth-order  Runge-Kutta  method  allows  times  steps  up  to  Courant  numbers  of  2  this 
ratio  represents  a  very  stringent  measure  of  the  Lagrange -Galerkin  time  step.  The  Courant 
number  on  unstructured  triangular  grids  is  best  obtained  in  an  edgewise  manner;  in  other 
words,  the  characteristic  length  used  is  the  sides  of  the  triangle  and  the  velocity  is  the  velocity 
vector  at  the  mid-point  of  the  particular  side,  thereby  yielding 


C  =  At 


/  u2  +  v2 
Ax2  +  Ay2 


It  should  also  be  noted  that  the  resolution  of  the  grid  is  given  in  terms  of  the  number  of 
quadrilaterals  in  the  grid.  For  instance,  a  40  x  40  grid  (TV  =  40)  actually  contains  80  x  80 
triangular  elements.  But  we  denote  the  resolution  of  the  grid  by  the  number  of  quadrilaterals 
as  these  yield  the  correct  element  spacing.  For  example,  for  a  domain  having  a  length  in  the 
v -direct ion  of  2  (as  in  cases  1  and  2),  the  TV=  40  grid  yields  an  element  spacing  of  Ax  =  2/40. 


6.1.  Problem  statement 

In  the  following  sections,  the  three  test  cases  used  to  measure  the  performance  of  the 
Lagrange -Galerkin  method  are  introduced.  They  are  a  solid  body  rotation  with  time -indepen¬ 
dent  velocity  field,  solid  body  rotation  with  time-dependent  velocity  field  and  a  westward 
traveling  soliton  wave. 

6.1.1.  Case  1:  solid  body  rotation.  The  initial  condition  for  this  test  case  is  given  by  the 
Gaussian  wave 


tp(  x,  y,  0)  =  exp 


-  [(  *  -  a-q)2  +  (y  -  Jo)2 
2A2 


having  the  far  boundary  conditions 
<p(x,  y,  t)  =  0,  V(x,  v)eT 
where 


1 


1 


4  =  o  >  Oo,  Vo)  =  (  -  2’  0  )>  <X  J)e[  -  1,  1] 
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The  velocity  field  is  constant  for  all  times  and  is  given  by 
u  =  +  y  and  v  =  —  x 


which  defines  a  clockwise  rotation  about  the  center  of  the  domain.  There  are  no  Coriolis 
effects  and  since  only  the  mass  is  allowed  to  vary,  this  problem  simplifies  to  the  passive 
advection  of  the  quantity  <p.  The  Gaussian  wave  rotates  along  this  circular  path  with  neither 
distortion  nor  dissipation.  The  exact  solution  is  given  by 


cp(x,  y,  t  )  =  exp 


-P)2  +  (J)2] 


where 


x  =  x  —  x0  cos  t  —  y0  sin  t  and  y  =  y  +  x0  sin  t  —  y0  cos  t 

Results  are  given  for  one  full  revolution  of  the  initial  wave.  The  period  for  one  revolution  of 
the  wave  is  2n,  which  means  that  one  revolution  corresponds  to  t  =  2n.  For  a  structured 
40  x  40  element  grid  (N  =  40)  ArRK  =  27t/80,  which  is  the  value  used  to  determine  the  time  step 
ratio  a  for  the  Lagrange-Galerkin  method. 

6.1.2.  Case  2:  solid  body  rotation  with  time -dependent  velocity.  This  test  case  was  introduced  by 
Chukapalli  [23],  The  only  difference  between  this  test  case  and  the  previous  one  is  that  the 
velocity  field  is  now  given  as 

u  =  +  v(  1  +  cos  t)  and  v  =  —  x(  1  +  cos  t )  (25) 

Therefore,  the  problem  has  the  same  period  of  revolution  as  in  case  1  and  it  admits  the  same 
type  of  solution  except  that  now  the  velocity  field  is  a  function  of  time.  Note  that  this  velocity 
field  is  uniformly  rotational  about  the  center  of  the  domain.  Therefore,  any  fluid  particle 
starting  out  at  a  specific  radial  distance  from  the  center  of  the  domain  will  always  remain  on 
the  circle  defined  by  the  initial  radius.  For  this  reason  it  is  best  to  obtain  the  analytic 
trajectories  by  mapping  from  Cartesian  to  polar  co-ordinates  using 

x  =  r  cos  9  and  y  =  r  sin  0  (26) 

Using  the  u  component  we  note  that 

d.v 

—  =  u  =  +  y{\  +  cos  t)  (27) 

dr 

and  substituting  (26)  into  the  left-hand  side  of  (27)  yields 
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d  dr  d  0 

—  (r  cos  6)  =  cos  6 - —  r  sin  0  — 
At  At  At 


But  note  that  from 


r  =  y.v2  +  v2 


we  can  show  that 


dr  x  d.r  y  Ay 

dr  r  dr  r  dr 


Substituting  the  u  and  v  components  of  velocity  from  Equation  (25)  yields 

d.r  x  y 

—  =  -y(l  +  cos  r)  —  -x(l  +cos  r)  =  0 

dr  r '  r 

Therefore,  Equation  (27)  becomes 

AO 

—  r  sin  0  —  =  +  y{\  +  cos  t ) 

dr 

and  canceling  like  terms  yields 

--  -(1+cosO 


(28) 


which  describes  the  trajectories  in  terms  of  the  angular  component  as  a  function  of  time  (we 
would  have  obtained  the  same  expression  if  we  had  used  the  v  component  instead).  Since  we 
know  that  the  exact  solution  is  of  the  form 


(p(x,  y,t  )  =  exp 


-  [(x  -  x0)2  +  (y  -  y0)2 


we  now  also  know  the  exact  solution  in  terms  of  0  because  the  Cartesian  co-ordinates  are  given 
in  terms  of  the  angular  co-ordinates  as 

x  =  r  cos  0(t)  and  y  =  y  sin  0(t) 

Integrating  Equation  (28)  we  get 

0(t)  =  0o  —  (r  +  sin  r  ) 

where 
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90  =  arctan  — 
x0 


This  test  problem  also  has  a  period  of  t  =  2n  but  the  velocity  field  initially  is  very  large  and 
slows  down  continuously  until  finally  reaching  zero  velocity  at  one  half  revolution  (t  =  n).  The 
velocity  field  then  begins  to  increase  at  this  point  and  by  virtue  of  the  behavior  of  the  cosine 
function,  the  velocity  field  is  symmetric  with  respect  to  the  time  t  =  n.  For  a  structured  40  x  40 
element  grid  (N  =  40),  A/RK  =  27r/160  for  this  case.  This  maximum  time  step  is  smaller  than  in 
case  1  owing  to  the  fact  that  the  velocity  field  is  twice  that  of  case  1  near  t  =  0  and  t  =  2n,  and 
therefore  a  time  step  half  as  large  is  required  for  stability. 

6.1.3.  Case  3:  Rossby  soliton  waves.  This  problem  describes  an  equatorially  trapped  Rossby 
soliton  wave  [27].  The  soliton  wave  starts  off  in  the  center  of  the  domain.  I  then  moves 
westward  along  the  equator  without  changing  shape.  The  exact  solution  is  given  by 

cp(x,  y,  t)  =  <p(0)  +  <pa> 

u(x,  y,  t )  =  m<0)  +  w(1) 

v(x,  y,  t )  =  v(0>  +  v(1> 


where  the  superscripts  (0)  and  (1)  denote  the  zeroth-  and  first-order  asymptotic  solutions  of  the 
shallow  water  equations  respectively.  They  are  given  by 


(p (0)  =  ij 


-9  +  6y 


-ri  2 


m(0)  =  xl  (2y)  e  ~yl/2 

' 


r(0:)  =  t] 


3  +  6y 


cpa)  =  c(1);/  —  (  —  5  +  2v2)  e  y  /2  +  ;/2Oa)(  v) 
16 


=  cah 7  —  (3  +  2 y2)  e  “■ ^/2  +  ;/2t7(1)(v) 
16 


(y) 

where  (c,  t)  =  A  sech2  Bq,  c  =  x  —  ct,  A  =  0.771R2,  B  =  0.394,  and  c  =  c(0)  +  c(1),  where  c(0)  = 
—  |  and  c(1)  =  —  0.39552.  The  variable  ;/  is  the  solution  to  the  equation 


Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Fluids  2000;  33:  789-832 


WEAK  LAGRANGE-GALERKIN  FINITE  ELEMENT  METHOD 


813 


fy  ,  Srj  d3// 

yT  +  an> l^  +  Pn-^3-0 


which  is  the  famous  Korteweg-de  Vries  equation,  which  yields  soliton  wave  solutions.  The 
shallow  water  equations  can  be  simplified  into  this  equation  using  the  method  of  multiple 
scales  presented  in  Reference  [28].  Finally,  the  remaining  terms  are  given  by 


oo 

r  x 

(Pn 

Um(y) 

=  e-F/2  £ 

u„ 

Vm(y) 

K  J 

n  —  0 

V„ 

V.  J 

where  Hn(y)  are  the  Hermite  polynomials  and  <pn,  un,  v„  are  the  Hermite  series  coefficients 
given  in  Reference  [27].  The  boundary  conditions  used  are 

n-u  =  0,  V(x,  v)eT 

where  n  is  the  outward  pointing  normal  vector  and  the  Coriolis  parameter  is  given  by 

f(y)  =  y 

where 

xe[  —  24,  +24]  and  ye[  —  8,  +8] 

(see  References  [1,2]  for  further  details  of  this  test  case).  The  equations  are  integrated  up  to  a 
non-dimensional  time  of  t  =  10.  For  a  structured  64  x  32  element  grid  AtRK  =  yg. 

6.2.  Results 

6.2.1.  Structured  grids 

6.2.1. 1.  Case  1.  This  case  was  run  using  structured  grids  as  shown  in  Figure  2.  This  plot  shows 
the  grid  and  contours  of  the  mass  for  a  40  x  40  element  grid  ( N  =  40).  Figure  3  shows  the  L2 
norm  as  a  function  of  the  number  of  elements  N  in  each  direction.  All  results  are  given  using 
80  time  steps  (A?  =  2n/80)  to  complete  one  full  revolution.  Therefore,  for  the  lowest  resolution 
grid  (AT  =20)  a  is  around  5,  while  for  the  highest  resolution  grid  (N  =  100)  a  is  close  to  3. 
Figure  3  shows  that  the  error  decreases  quadratically  with  N. 

Figure  4  shows  the  effects  of  a  on  the  solution  of  an  /V  =  40  grid  for  the  following  three 
trajectory  calculation  methods:  Runge-Kutta  (RK),  composite  mid-point  rule  (CM),  and 
exact  trajectories  (EX).  Specifically,  Figure  4(a)  shows  the  effects  on  the  L2  norm  of  the  mass 
while  Figure  4(b)  shows  the  trajectory  norm.  From  Figure  4(a)  we  see  that  the  RK  method  is 
more  accurate  than  the  CM  method.  Clearly,  using  exact  trajectories  yields  far  more  accurate 
results  but  they  are  hardly  ever  known  a  priori  except  in  transport  models  where  the  velocity 
fields  are  known  at  all  times.  Figure  4(b)  shows  that  the  trajectory  error  increases  linearly  with 
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Figure  2.  The  IV  =  40  structured  grid  and  mass  contours  for  cases  1  and  2  at  t  =  2n. 

the  time  step.  However,  Figure  4(a)  shows  that  there  is  a  trade-off  between  time  step  size  and 
the  number  of  interpolations.  In  other  words,  for  small  time  steps  very  accurate  trajectories  are 
obtained  but  many  interpolations  are  required,  while  large  time  steps  yield  less  accurate 
trajectories  but  require  fewer  interpolations. 

This  test  case  was  also  used  to  compare  the  piecewise  exact  and  numerical  integration 
approaches  for  obtaining  the  integrals  at  the  feet  of  the  characteristics  (the  integrals  involving 
the  Lagrangian  elements,  i.e.,  those  at  time  tn).  Numerical  experiments  were  performed  on  the 
N  =  40  grid  and  the  exact  and  numerical  approaches  yielded  error  norms  within  0. 1  percent  of 
each  other  with  similar  conservation  measures  M  and  E;  however,  the  exact  integration 
approach  took  twice  as  long  to  complete  one  full  revolution.  The  inefficiency  of  the  algorithm 
is  in  the  decomposition  of  the  Lagrangian  element  into  the  intersection  regions  between  the 
Lagrangian  and  Eulerian  elements,  which  are  the  regions  having  the  basis  functions  given  in 
Equation  (23).  An  advancing  front  mesh  generator  and  a  point  insertion  Delaunay  triangula¬ 
tion  were  both  compared  and  although  the  Delaunay  method  was  much  faster  it  still  could  not 
compete  with  the  numerical  integration  approach.  On  these  grounds,  the  numerical  integration 
approach  is  used  exclusively  for  the  remainder  of  the  paper. 

6.2. 1.2.  Case  2.  This  test  case  is  only  used  to  test  the  accuracy  of  the  trajectory  calculation 
methods.  In  case  1  we  were  able  to  test  the  trajectory  calculations  for  a  steady  velocity  field. 
Case  2  is  used  to  test  the  trajectory  calculations  for  a  transient  velocity  field.  A  comparison  of 
the  CM  and  RK  schemes  is  shown  in  Figure  5  for  a  structured  N  =  40  grid.  Figure  5(a)  shows 
the  effects  of  varying  the  time  step  ratio  a  on  the  L2  norm  of  the  mass  while  Figure  5(b)  shows 
the  effects  on  the  trajectory  norm.  It  is  evident  from  this  figure  that  the  RK  method  is  superior 
to  the  CM  method  but  now  only  slightly  so.  It  is  quite  surprising  that  the  difference  between 
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Figure  3.  The  mass  L2  norm  as  a  function  of  the  number  of  structured  elements  N  for  case  1  at  t  =  2k. 

the  RK  and  CM  methods  is  so  slight  especially  since  Figure  5(b)  shows  that  the  RK  method 
yields  far  more  accurate  trajectory  norms  than  the  CM  method.  Therefore,  for  such  a  strongly 
varying  velocity  field  as  used  in  this  test  case,  most  of  the  error  to  the  mass  comes  from  the 
extrapolation  of  the  velocity  field  and  not  from  the  trajectory  calculation.  Figure  5(a)  shows 
that  if  the  exact  trajectories,  or  values  close  to  them,  are  used  then  a  very  good  error  norm  for 
the  mass  can  be  obtained.  Clearly,  there  is  room  for  improvement  and  although  we  may  not 
expect  to  obtain  the  exact  trajectories,  we  should  be  able  to  come  quite  close. 

6.2. 1.3.  Case  3.  Figure  6  shows  the  grid  and  contours  for  a  structured  64  x  32  (N  =  64)  element 
grid  using  100  iterations  to  get  to  t  =  10.  Figure  6  shows  the  (a)  grid,  (b)  mass,  (c)  u 
momentum,  and  (d)  v  momentum.  Since  the  flow  is  unidirectional  along  the  x  co-ordinate,  the 
v  momentum  is  actually  quite  small  and  very  close  to  zero.  Qualitatively,  we  would  like  the 
solution  to  be  symmetric  with  respect  to  the  y-axis.  Figure  7  shows  the  L2  norm  for  the  three 
primitive  variables  as  a  function  of  N.  For  this  test  case,  N  refers  to  the  number  of  elements 
in  the  x -direction.  The  number  of  elements  in  the  y -direction  is  Nj 2.  Therefore,  for  an  N  =  80 
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Figure  5.  The  effect  of  the  time  step  ratio  a  on  the  (a)  mass  L2  norm  and  (b)  trajectory  norm  T  for  the 
structured  N  =  40  grid  for  case  2  at  t  =  2n  using  trajectories  obtained  by:  the  RK  method,  the  CM  rule, 

and  EX. 
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Figure  6.  The  structured  iV=64  (a)  grid  and  contours  for  the  (b)  mass,  (c)  u  momentum,  and  (d)  v 

momentum  for  case  3  at  t=  10. 
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Figure  7.  The  L2  norms  as  a  function  of  the  number  of  structured  elements  N  for  case  3  at  t  =  10. 

grid,  Nx  =  80  and  Ny  =  40,  and  so  on.  In  Figure  7,  the  mass  and  u  momentum  are  converging 
at  a  quadratic  rate,  while  the  v  momentum  is  converging  at  a  slightly  lower  rate. 

Figure  8  shows  the  effects  of  the  time  step  ratio  a  on  the  L2  norm  using  an  N  =  64  grid.  The 
error  norms  remain  the  same  for  almost  all  values  of  er,  except  around  5,  where  the  solution 
accuracy  begins  to  diminish.  Therefore,  as  was  seen  in  cases  1  and  2,  the  Lagrange- Galerkin 
method  has  a  clear  optimal  time  step.  It  would  be  very  useful  to  be  able  to  determine  what  this 
optimal  time  step  should  be;  however,  there  is  definitely  a  maximum  time  step  that  should  not 
be  exceeded.  Increasing  the  time  step  beyond  this  maximum  allowable  time  step  greatly 
diminishes  the  numerical  accuracy  of  the  solution. 

6.2.2.  Unstructured  grids 

6.2.2. 1 .  Case  1.  This  case  was  run  using  unstructured  grids  like  those  shown  in  Figure  9.  This 
figure  illustrates  the  grid  and  contours  of  the  mass  for  a  comparable  grid  to  the  structured 
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Figure  8.  The  effect  of  the  step  ratio  a  on  the  L2  norms  for  the  structured  N  =  64  grid  for  case  3  at 

(=10. 

#  =  40  grid.  Figure  10  shows  the  L2  norm  as  a  function  of  N.  All  results  are  obtained  using 
80  time  steps  to  complete  one  full  revolution.  Therefore,  for  the  lowest  resolution  grid  (N  =  10) 
a  is  also  around  and  for  the  highest  resolution  grid  (N  =  100)  a  is  around  3.  Figure  10  shows 
that  the  error  decreases  quadratically  with  N  even  for  these  unstructured  grids.  In  fact,  by 
comparing  Figures  3  and  10  we  can  see  that  both  the  structured  and  unstructured  grids  yield 
about  the  same  order  of  accuracy.  This  case  shows  that  the  weak  Lagrange -Galerkin  method 
works  equally  well  on  both  structured  and  unstructured  grids. 

6.2.22.  Case  3.  The  results  presented  for  this  case  are  obtained  using  100  iterations  to  reach 
the  non-dimensional  time  of  t  =  10;  however,  the  grids  are  now  unstructured.  These  grids  are 
not  randomly  unstructured  but  rather  the  very  flexible  mesh  generation  scheme  presented  here 
is  used  to  construct  a  priori  grids.  A  priori  grids  refer  to  the  construction  of  a  grid  by 
anticipating  the  behavior  of  the  solution.  In  this  case,  since  the  soliton  wave  moves  along  a 


Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Fluids  2000;  33:  789-832 


WEAK  LAGRANGE-GALERKIN  FINITE  ELEMENT  METHOD 


821 


Figure  9.  The  N  =  40  unstructured  grid  and  mass  contours  for  case  1  at  t  =  2n. 

very  constrained  region  along  the  x  -direction,  it  makes  sense  to  introduce  many  elements 
where  the  wave  is  destined  to  pass  and  few  elements  where  the  wave  will  not.  Figure  1 1  shows 
a  sample  of  an  a  priori  grid  comparable  with  the  structured  N  =  64  grid.  Figure  1 1  shows  the 
(a)  grid,  (b)  mass,  (c)  u  momentum,  and  (d)  v  momentum.  Figure  12  shows  the  L2  norm  as  a 
function  of  N.  The  plot  shows  a  fast  rate  of  convergence  especially  for  the  mass  (<p)  and  u 
momentum  («).  This  test  case  shows  qualitatively  and  quantitatively  that  the  weak  Lagrange - 
Galerkin  method  yields  very  good  solutions  to  the  full  non-linear  shallow  water  equations  on 
unstructured  grids.  In  fact,  the  unstructured  grids  achieved  the  same  order  of  accuracy  as  the 
structured  grids  as  is  evidenced  by  comparing  Figures  7  and  12. 

6.2.3.  Adaptive  grids 

6.2.3. 1.  Case  1.  As  a  final  test,  this  case  was  run  using  an  adaptive  unstructured  mesh 
generation  strategy.  The  initial  grid  is  comparable  with  the  structured  N  =  40  grid  and  is 
allowed  to  adapt  only  after  each  quarter  revolution;  however,  the  number  of  points  allowed 
are  restricted  to  give  a  grid  comparable  with  the  N  =  40  structured  grid.  As  in  the  previous 
N  =  40  cases,  80  iterations  were  taken  to  reach  one  full  revolution  of  the  Gaussian  wave. 
Figure  13  shows  the  grid  and  contours  at  (a)  t  =  n/2,  (b)  t  =  n,  (c)  t  =  3n/2,  and  (d)  t  =  2n. 
Note  that  the  mesh  generator  captured  the  Gaussian  wave  extremely  well,  which  then  assists 
the  Lagrange-Galerkin  method  in  resolving  the  wave  quite  sharply.  In  fact,  we  get  a  50 
percent  increase  in  accuracy  by  using  this  strategy  over  the  structured  and  unstructured  grids. 

6.2. 3.2.  Case  3.  Figure  14  shows  the  adaptive  grids  and  the  mass  contours  for  the  soliton  wave 
at  times  (a)  t  =  2.5,  (b)  t  =  5,  (c)  t  =  7.5,  and  (d)  t  =  10.  Once  again  100  iterations  were  used 
to  integrate  the  solution.  The  adaptive  grids  used  are  comparable  with  the  N  =  64  structured 
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Figure  10.  The  mass  L2  norm  as  a  function  of  the  number  of  unstructured  elements  N  for  case  1  at 

t  =  2n. 


grid  and  were  allowed  to  adapt  only  at  the  time  given  above.  The  solution  has  been 
successfully  tracked  by  both  the  Lagrange -Galerkin  method  and  the  adaptive  mesh  generator. 
Figure  15  shows  the  (a)  grid,  (b)  mass,  (c)  u  momentum,  and  (d)  v  momentum  at  time  t  =  10. 
Figure  15(a)  clearly  shows  two  distinct  wave  structures,  one  on  top  of  the  other.  Figure  15(b) 
shows  the  mass  solution  and  there  are  clearly  two  distinct  wave  structures  represented  here  as 
well.  Finally,  Figure  15(c)  and  (d)  show  the  u  and  v  momentum  respectively.  The  u  momentum 
looks  very  much  symmetric  with  respect  to  the  y-axis  and  the  v  momentum  is  more  or  less 
symmetric,  however,  these  values  are  extremely  small,  as  the  flow  is  essentially  unidirectional 
along  the  x-axis.  This  test  case  truly  shows  the  power  and  flexibility  of  using  the  weak 
Lagrange -Galerkin  method  in  conjunction  with  adaptive  unstructured  grids.  In  fact,  the 
accuracy  of  this  solution  over  the  structured  and  unstructured  grids  is  around  25  percent  more 
accurate. 
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Figure  11.  The  unstructured  N  =  64  (a)  grid  and  contours  for  the  (b)  mass,  (c)  u  momentum,  and  (d)  v 

momentum  for  case  3  at  t=  10. 
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Figure  12.  The  L2  norms  as  a  function  of  the  number  of  unstructured  elements  N  for  case  3  at  t  =  10. 

The  reason  why  Lagrangian  methods  are  so  desirable  for  use  with  adaptive  mesh  generators 
is  that  the  mesh  generator  constructs  very  small  elements  in  certain  regions,  which  then  alters 
the  maximum  local  Courant  number.  If  Eulerian  methods  were  used,  the  global  time  step 
would  have  to  be  changed  in  order  to  satisfy  the  CFL  condition  at  these  small  elements.  This 
is  true  not  just  for  explicit  methods  (due  to  instabilities)  but  for  implicit  methods  (due  to 
inaccuracies)  as  well  since  the  goal  of  using  adaptive  grids  is  to  achieve  high  accuracy  as  well 
as  stability;  Lagrangian  methods  offer  both. 

The  other  attraction  of  using  adaptive  grids  is  that  we  do  not  have  to  concern  ourselves  with 
trying  to  predict  how  the  solution  will  behave.  Instead,  we  feed  the  solution  algorithm  a  crude 
mesh  and  then  let  the  Lagrange -Galerkin  method  and  the  adaptive  mesh  generator  capture 
the  interesting  phenomena  of  a  given  problem.  However,  the  initial  grid  must  be  sufficiently 
fine  such  that  it  can  capture  the  small-scale  features  of  the  flow.  In  other  words,  assume  we 
have  a  high  resolution  adapted  grid  at  time  t".  In  order  to  maintain  the  same  level  of  accuracy 
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Figure  14.  The  adaptive  unstructured  N  =  64  grid  and  mass  contours  at  (a)  t  =  2.5,  (b)  t  =  5,  (c)  t  =  7.5 

and  (d)  t=  10  for  case  3. 


at  time  tn  +  1  we  require  a  fine  mesh  where  the  wave  will  move  to  because  if  we  have  a  crude 
mesh  there,  then  no  matter  how  good  the  Lagrange -Galerkin  method  is  it  will  not  be  able  to 
maintain  the  same  level  of  accuracy.  For  situations  where  the  velocity  field  is  known  at  all 
times,  such  as  in  air  pollution  transport,  we  can  use  adaptive  grids  at  both  ends  of  the 
characteristics  (times  t"  and  tn+ ')  in  order  to  get  very  precise  solutions.  This  idea  can  also  be 
used  to  track  certain  tracers,  fumes,  or  plumes. 

6.3.  Computational  cost 

Figure  16  summarizes  the  computational  cost  of  using  the  Lagrange-Galerkin  method  for  the 
TV  =64  structured  grid  for  case  3.  These  results  were  obtained  on  a  DEC  Alpha  8400  (EV5 
CPU  chip)  running  only  one  processor  at  a  clock  speed  of  300  MHz.  The  computational  cost 
is  divided  into  three  categories:  solve,  search,  and  fem.  Solve  represents  the  percent  of  the  total 
time  to  invert  the  mass  and  momentum  matrices.  This  includes  the  LU  decompositions  and  the 
back  substitution  of  these  matrices.  Search  involves  the  percent  of  time  required  for  the 
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Figure  15.  The  adaptive  unstructured  A^=64  (a)  grid  and  contours  for  the  (b)  mass,  (c)  u  momentum, 

and  (d)  v  momentum  for  case  3  at  t  =  10. 
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departure  points  calculations,  quadtree  searches,  and  testing  for  inclusions.  Fern  denotes  the 
percent  involving  any  operation  typically  associated  with  the  finite  element  method.  This 
includes  the  construction  of  the  left-  and  right-hand  side  element  equations  as  well  as  the 
summation  of  these  element  contributions,  which  result  in  the  global  system  of  equations.  The 
percentage  breakdown  shown  in  Figure  16  is  based  on  a  1170-s  computation. 

The  interesting  thing  about  Figure  16  is  that  neither  the  searching  algorithms  nor  the  finite 
element  operations  take  any  time  whatsoever  relative  to  the  amount  of  time  required  to  invert 
the  global  matrices.  What  this  means  is  that  if  the  algorithm  must  be  made  more  efficient,  it 
will  only  be  accomplished  by  either  constructing  a  more  efficient  matrix  solver  or  by  picking 
a  better  one. 

Many  other  possibilities  for  increasing  the  efficiency  of  the  current  algorithm  can  be 
explored:  multi-grid-type  methods  have  proven  to  be  the  most  efficient  for  solving  linear 
equations,  the  LU  solver  can  be  streamlined  by  using  a  node  renumbering  scheme,  which 
would  then  allow  the  storage  of  the  matrix  in  a  tightly  banded  form,  and  finally  a  parallel 
implementation  of  the  existing  algorithms  would  remove  any  of  the  current  bottlenecks. 
Currently  we  are  exploring  multi-grid  methods  and  parallel  implementations  using  the  message 
passing  interface  (MPI). 


Figure  16.  The  percentage  breakdown  of  the  computational  costs  incurred  by  the  various  operations  of 
the  Lagrange -Galerkin  method  for  the  structured  N  =  64  grid  for  case  3  at  t  =  10.  The  total  CPU  time 

consisted  of  1170  s. 
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7.  CONCLUSIONS 

The  weak  Lagrange-Galerkin  finite  element  method  on  adaptive  unstructured  grids  is 
presented.  In  order  to  apply  the  weak  Lagrange -Galerkin  method  to  the  shallow  water 
equations  they  must  first  be  written  in  conservation  form.  The  combination  of  the  conservative 
weak  Lagrange-Galerkin  method  along  with  the  governing  equations  having  been  written  in 
conservation  form  results  in  a  highly  conservative  scheme.  In  fact,  the  mass  and  energy  were 
exactly  conserved  (up  to  99.8  percent)  for  all  the  grids  and  for  all  of  the  test  cases  regardless 
of  the  size  of  the  time  step.  The  weak  Lagrange-Galerkin  method  produces  a  very  efficient 
scheme  for  solving  the  shallow  water  equations  that  are  not  just  applicable  to  the  two-dimen¬ 
sional  planar  shallow  water  equations  but  are  also  directly  extendible  to  the  two-dimensional 
spherical  and  three-dimensional  Cartesian  shallow  water  equations.  This  generalization  is 
possible  through  the  use  of  the  natural  co-ordinates  as  the  finite  element  basis  functions. 

Although  the  shallow  water  equations  are  a  coupled  system  of  non-linear  first-order  PDEs, 
after  discretizing  by  the  weak  Lagrange-Galerkin  method  they  yield  a  decoupled  system  of 
linear  first-order  equations.  The  decoupling  occurs  between  the  mass  and  the  momentum 
equations;  the  n  momentum  equations  (for  R"  space)  remain  weakly  coupled  through  the 
Coriolis  term.  In  this  paper  we  have  described  in  detail  the  application  of  the  weak 
Lagrange-Galerkin  method  to  the  shallow  water  equations.  In  additions,  we  have  shown  how 
to  integrate  all  of  the  equations  exactly,  not  just  those  at  the  arrival  points  (the  Eulerian 
elements  at  time  t"  +  1)  but  also  the  departure  points  (the  Lagrangian  elements  at  time  t"). 
However,  it  was  found  that  the  numerical  integration  approach  yielded  results  equal  or  close 
to  the  exact  integration  approach  while  costing  a  factor  of  two  less.  But  since  this  is  such  and 
elegant  method,  further  work  to  increase  its  efficiency  needs  to  be  explored. 

Perhaps  one  of  the  most  important  aspects  of  a  Lagrangian  method  is  the  correct  capturing 
of  the  fluid  particle  trajectories.  In  this  paper,  we  have  introduced  two  new  methods  for 
calculating  these  trajectories:  the  Runge-Kutta  (RK)  and  the  composite  mid-point  rule  (CM). 
Both  methods  are  shown  to  be  extremely  accurate  not  just  for  steady  state  velocity  fields  but 
for  transient  ones  as  well.  The  RK  method  performed  better  than  the  CM  method  but  further 
improvements  can  always  be  made. 

Linear  triangular  elements  are  used  in  this  study,  which  give  three  advantages  over 
quadrilaterals:  exact  or  numerical  integration  can  be  used  to  integrate  all  of  the  finite  element 
integrals;  regardless  of  shear  or  rotational  aspects  of  the  flow,  triangles  will  not  become 
non-convex  polygons  whereas  this  danger  exists  and  is  very  probable  with  quadrilaterals;  and 
finally  adaptive  unstructured  grids  can  be  used.  The  adaptive  unstructured  grid  generator 
presented  is  an  advancing  front  method,  which  uses  the  idea  of  the  weak  Lagrange-Galerkin 
method  to  project  the  values  at  the  grid  points  from  the  old  mesh  onto  the  newly  adapted 
mesh.  This  idea  results  in  a  difference  in  conservation  of  mass  and  energy  from  one  mesh  to 
the  other  on  the  order  of  1  x  10  _ 8  percent. 

Three  test  cases  were  run  in  order  to  validate  the  algorithms:  a  solid  body  rotation  of  a 
Gaussian  with  a  steady  state  velocity  field,  a  solid  body  rotation  of  a  Gaussian  with  a  transient 
velocity  field,  and  a  westward  traveling  soliton  wave.  The  first  two  test  cases  are  given 
primarily  to  show  that  the  convergence  rate  of  the  weak  Lagrange-Galerkin  method  is 
quadratic  (second-order)  and  to  compare  the  two  methods  for  calculating  the  fluid  particle 
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trajectories.  The  soliton  case  tests  the  ability  of  the  scheme  to  handle  a  non-linear  phenomenon 
and  it  performed  extremely  well,  not  only  for  the  structured  and  unstructured  grids  but  also 
for  the  adaptive  unstructured  grids. 

The  results  obtained  for  the  three  problems  demonstrate  that  the  weak  Lagrange-Galerkin 
finite  element  method  with  adaptive  unstructured  triangles  offer  a  viable  method  for  solving 
the  shallow  water  equations  on  the  plane.  This  paper  describes  the  approach  in  detail  while 
future  papers  will  focus  on  the  application  of  this  method  for  various  types  of  problem.  The 
development  of  a  three-dimensional  weak  Lagrange-Galerkin  shallow  water  model  is  cur¬ 
rently  underway. 
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APPENDIX  A 


In  order  to  show  that  the  basis  functions  vanish  along  the  characteristics,  let  us  assume  for 
simplicity  that  we  can  write  the  basis  functions  as  a  linear  function  of  time.  Therefore,  between 
times  ?"  and  f'!+1  we  can  write  the  basis  functions  as 


<A,(x,  t) 


lA?+1(x)  + 


<^(x) 


(29) 


recalling  that  the  basis  functions  at  tn+1  (i.e.,  at  the  Eulerian  elements,  see  Figure  1)  are 


^+1(x) 


¥  +  bty  +  Cj 

det 


(30) 


where 

det  AUJc  =  (xj  -  x,.)  x  (x*  -  x,) 


and 


at  =  V]  -  Tfo  b,  =  xk  -  Xj,  c,  =  Xjyk  -  xkyj 

Assuming  a  constant  velocity  field,  from  the  trajectories  we  can  write 
x"  =  Xj  —  a  At  and  y”  =  yt  —  v  At 
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where 


<"  +  1  =  x,-  and  y" + 1  =  y, 


Plugging  in  the  x"  values  into  Equation  (30)  and  simplifying  gives 

atu  +  b,v 


^(x)  =  ^  +  1(x)  +  At 


det  A,- 


(31) 


i,j,k 


Substituting  Equations  (30)  and  (31)  into  Equation  (29)  and  grouping  like  terms  yields 

'  a, u  +  b.v ' 


^.(x,  0  =  ^"+1(x)  +  (f"+1-0 

Differentiating  yields 


del  A  fjjc 


and 


f  atu  +  btv 

8t 

<f 

<D 

_3 

dh_ 

f  ai 

8x 

Idct  A  iJJc 

di//i 

(  b‘ 

Sy 

\det  A  ijk 

Substituting  into  the  characteristic  equation 

#«■  #«,  #<■  n 

of  ox  op 

we  can  see  that  the  basis  functions  do  indeed  vanish  along  the  characteristics.  A  similar 
analysis  can  be  carried  out  for  the  general  case  where  u  and  v  are  no  longer  constant;  however, 
it  is  considerably  more  complicated. 
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